function [va_t, va_v] = velocityanalyse(shotw, th, dt, velrange, nvel, nt0, win, smoothL)

% stack
t0range = [0, (size(shotw, 1)-1)*dt];
xrecs = arrayfun(@(x) x.offset, th);
va = zeros(nt0, nvel);
tl = (0:size(shotw, 1)-1)*dt;
infostep = round((0:0.2:1)*(nvel-1))+1;
seismap = interp1([0, 1, 32, 64], [1,1,1; 0,0,0;1,1,1;1,0,0], linspace(-1, 1, 65).^3*32+32);
tic;
for j = 1:nvel
    if any(infostep == j)
        fprintf('stack progress: %.0f%%\n', 100*(j-1)/(nvel-1));
    end
    vel = velrange(1) + (j-1)/(nvel-1)*(velrange(2)-velrange(1));
    shot_s = nmor(shotw, (0:size(shotw, 1)-1)*dt, xrecs, vel, 1);
%   shot_s = normalmoveout(shotw, xrecs, dt, 0, vel);
    parfor i = 1:nt0
        t0 = t0range(1)+(i-1)/(nt0-1)*(t0range(2)-t0range(1));
        cutwin = shot_s(tl>t0-win & tl < t0+win, :);
%         cutwin = stackvel(shotw, xrecs, dt, t0-win:dt:t0+win, 0, vel);
        va(i, j) = sum(sum(cutwin, 2).^2, 'omitnan')/sum(cutwin(:).^2, 'omitnan');
    end
end
toc;
% smooth
sm_x = (-smoothL:smoothL)/smoothL;
ker = exp(-sqrt(sm_x.^2+sm_x'.^2)*2);
va = conv2(va, ker);
va = va/max(va(:));
% va = va.^2;
% pick
va_v = [];
va_t = [];
figure('Position', [10 100 1887 838], 'Units', 'normalized');
ax1 = axes('Position', [0.05, 0.05, 0.4, 0.9]);
imagesc(ax1, velrange, t0range, va);hold on;shading interp;colormap(parula);
ax2 = axes('Position', [0.45, 0.05, 0.2, 0.9]);
plotseiswave(shotw, [th.offset], dt, 'trace', 'line', ax2);
ax3 = axes('Position', [0.65, 0.05, 0.1, 0.9]);
plotseiswave(shotw, [th.offset], dt, 'trace', 'image', ax3);
colormap(ax3, seismap);
ax4 = axes('Position', [0.75, 0.05, 0.1, 0.9], 'Box', 'on');
ax5 = axes('Position', [0.85, 0.05, 0.1, 0.9], 'Box', 'on');
t = 0:360;
r = 0.01;
rx = r*(velrange(2)-velrange(1));
ry = r*(t0range(2)-t0range(1));
while true
    % pick
    [x, y, b] = ginput(1);
    if b == 2
        break
    elseif b == 3
        dist = sqrt((x-va_v).^2/rx^2+(y-va_t).^2/ry^2);
        idx = find(dist>1);
        if isempty(idx)
            va_v = [];
            va_t = [];
        else
            va_v = va_v(idx);
            va_t = va_t(idx);
        end
    else
        va_v = [va_v, x];
        va_t = [va_t, y];
    end
    % sort by time
    [~, idx] = sort(va_t, 'ascend');
    va_t = va_t(idx);
    va_v = va_v(idx);
    % plot 1
    for ax = [ax1, ax2, ax3, ax4, ax5]
        cla(ax);
    end
    imagesc(ax1, velrange, t0range, va);hold on;shading interp;colormap(parula);
    line(ax1, va_v, va_t, 'Color', 'r', 'LineWidth', 1);
    if ~isempty(va_t)
        for i = 1:length(va_t)
            line(ax1, va_v(i)+cosd(t)*rx, va_t(i)+sind(t)*ry, 'Color', 'r', 'LineWidth', 1);
        end
    end
    % plot 2 3
    if length(va_v)<1
        plotseiswave(shotw, [th.offset], dt, 'trace', 'line', ax2);
        plotseiswave(shotw, [th.offset], dt, 'trace', 'image', ax3);
    else
        shot_s = nmor(shotw, (0:size(shotw, 1)-1)*dt, xrecs, fillvel(va_t, va_v, tl), 1);
        plotseiswave(shot_s, [th.offset], dt, 'trace', 'line', ax2);
        plotseiswave(shot_s, [th.offset], dt, 'trace', 'image', ax3);
    end
    colormap(ax3, seismap);
    if ~isempty(va_t)
        % plot 4
        line(ax4, fillvel(va_t, va_v, tl), tl, 'Color', 'k', 'LineWidth', 2);
        % plot 5
        [tint, vint] = inverseDixFormula(va_t, va_v);
        line(ax5, fillvel(tint, vint, tl), tl, 'Color', 'k', 'LineWidth', 2);
    end
    for ax = [ax2, ax3, ax4, ax5]
        ax.Box='on';
        ax.YDir='reverse';
        ax.YTickLabel={};
    end
    for ax = [ax2, ax3]
        ax.XTickLabel={};
    end
end
end

function y = stackvel(data, x, dt, tcut, vel_t, vel_v)
tline = (0:size(data, 1)-1)'*dt;
y = zeros(length(tcut), size(data, 2));
t0 = mean(tcut);
vcut = fillvel(vel_t, vel_v, t0);
for j = 1:size(y, 2)
    t = sqrt(t0.^2+(x(j)./vcut).^2);
    y(:, j) = interp1(tline, data(:, j), t+tcut-mean(t0), 'linear');
end
end